//============================================================================
// Name        : csiena
// Author      : Josh Lospinoso
// Copyright   : MIT License
// Description : Austere, limited implementation of the R package RSiena in
//				 C++. Used for fast maximum likelihood and method of moments
//				 estimation of stochastic actor oriented models (SAOMs)
//				 without any frills or user experience niceties.
//============================================================================

#include "crimp.h"

namespace csiena {

ChainCrimp::~ChainCrimp() {
	#if GAMMA_LENGTH > 0
	delete(this->score);
	delete(this->information);
	delete(this->statistic);
	#endif
	delete(this->observedDigraph);
}

ChainCrimp::ChainCrimp(int actors) {
	this->actors = actors;
	this->previous=NULL;
	this->next=NULL;
	this->logProbability=0;
	#if GAMMA_LENGTH > 0
		this->score = new double[GAMMA_LENGTH];
		this->information  = new double[GAMMA_LENGTH*GAMMA_LENGTH];
		this->statistic = new double[GAMMA_LENGTH*2];
	#endif
	this->observedDigraph = NULL;
	for (int k=0; k<GAMMA_LENGTH; k++) {
		this->setScore(k,0);
		for (int l=0; l<GAMMA_LENGTH; l++) {
			this->setInformation(k,l,0);
		}
	}
	for (int k=0; k<GAMMA_LENGTH*2; k++) {
		this->setStatistic(k,0);
	}
}

void ChainCrimp::evaluate(Digraph* x, Evaluation* fIn) {
	CrimpEvaluation* f = (CrimpEvaluation*)fIn;
	f->evaluate(this->observedDigraph, x);
	this->setLogProbability(f->getLogLikelihood());
	for (int k=0; k<GAMMA_LENGTH; k++) {
		this->setScore(k,f->getScore(k));
		for (int l=0; l<GAMMA_LENGTH; l++) {
			this->setInformation(k,l, f->getInformation(k,l));
		}
	}
	for (int k=0; k<GAMMA_LENGTH*2; k++) {
		this->setStatistic(k,f->getStatistics(k));
	}
}

double ChainCrimp::evaluateProposal(Digraph* x, Evaluation* fIn) {
	CrimpEvaluation* f = (CrimpEvaluation*)fIn;
	f->evaluate(this->observedDigraph, x);
	return f->getLogLikelihood();
}

void ChainCrimp::simulate(Digraph* trueDigraph, Digraph* observedDigraph, CrimpEvaluation* f, MersenneTwister* r) {
	if (GAMMA_LENGTH == 0) {
		for (int i=0; i<this->actors; i++) {
			for (int j=0; j<this->actors; j++) {
				observedDigraph->set(i,j, trueDigraph->get(i,j));
			}
		}
	} else if (GAMMA_LENGTH == 1) {
		for (int i=0; i<this->actors; i++) {
			for (int j=0; j<this->actors; j++) {
				if (r->next() > f->getGamma(0)) {
					observedDigraph->set(i,j, trueDigraph->get(i,j));
				} else {
					observedDigraph->set(i,j, !trueDigraph->get(i,j));
				}
			}
		}
	} else if (GAMMA_LENGTH == 2) {
		for (int i=0; i<this->actors; i++) {
			for (int j=0; j<this->actors; j++) {
				double errorRate = f->getGamma(1); /* False Positive Rate */
				if (trueDigraph->get(i,j)) {
					errorRate = f->getGamma(0); /* False Negative Rate */
				}
				if (r->next() > errorRate) {
					observedDigraph->set(i,j, trueDigraph->get(i,j));
				} else {
					observedDigraph->set(i,j, !trueDigraph->get(i,j));
				}
			}
		}
	} else {
		std::cerr << "Gamma = " << GAMMA_LENGTH << " not specified for simulation!" << std::endl;
	}
}

// Utility methods
void ChainCrimp::print(int index) {
	printf(" (%3u) CRIMP %+6.3f ", index, this->getLogProbability());
	printf(" |");
	for (int k=0; k<GAMMA_LENGTH; k++) {
		printf(" %+3.2f", this->getScore(k));
	}
	printf(" |");
	for (int k=0; k<GAMMA_LENGTH; k++) {
		printf(" %+3.2f", this->getInformation(k,k));
	}
	printf(" |");
	for (int k=0; k<GAMMA_LENGTH * 2; k++) {
		printf(" %+3.2f", this->getStatistic(k) );
	}
	printf("\n");
}

bool ChainCrimp::isCrimp() {
	return true;
}

// Accessor Methods
double ChainCrimp::getStatistic( int statisticIndex) {
	if (statisticIndex < 2*GAMMA_LENGTH) {
		return this->statistic[ statisticIndex];
	} else {
		std::cout << "ChainCrimp::getStatistic: index out of bounds." << statisticIndex << std::endl;
		return 0.0;
	}
}

void ChainCrimp::setStatistic( int statisticIndex, double value) {
	if (statisticIndex < 2*GAMMA_LENGTH) {
		this->statistic[statisticIndex] = value;
	} else {
		std::cout << "ChainCrimp::setStatistic: index out of bounds." << statisticIndex << std::endl;
	}
}

double ChainCrimp::getScore(int gammaIndex) {
	if (gammaIndex < GAMMA_LENGTH) {
		return this->score[gammaIndex];
	} else {
		std::cout << "ChainCrimp::getScore: index out of bounds." << gammaIndex << std::endl;
		return 0.0;
	}
}

void ChainCrimp::setScore(int gammaIndex, double value) {
	if (gammaIndex < GAMMA_LENGTH) {
		 this->score[gammaIndex] = value;
	} else {
		std::cout << "ChainCrimp::setScore: index out of bounds." << gammaIndex << std::endl;
	}
}

double ChainCrimp::getInformation(int gammaIndex1, int gammaIndex2) {
	if ((gammaIndex1 < GAMMA_LENGTH) & (gammaIndex2 < GAMMA_LENGTH)) {
		return this->information[gammaIndex1*GAMMA_LENGTH + gammaIndex2];
	} else {
		std::cout << "ChainCrimp::getInformation: index out of bounds " << gammaIndex1 << "," << gammaIndex2 << std::endl;
		return 0.0;
	}
}

void ChainCrimp::setInformation(int gammaIndex1, int gammaIndex2, double value) {
	if ((gammaIndex1 < GAMMA_LENGTH) & (gammaIndex2 < GAMMA_LENGTH)) {
		 this->information[gammaIndex1*GAMMA_LENGTH + gammaIndex2] = value;
	} else {
		std::cout << "ChainCrimp::setInformation: index out of bounds " << gammaIndex1 << "," << gammaIndex2 << std::endl;
	}
}

double ChainCrimp::getLogProbability() {
	return this->logProbability;
}

void ChainCrimp::setLogProbability(double value) {
	this->logProbability = value;
}

}
